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Non-equilibrium molecular dynamics simulations, of crucial importance in sliding friction, are 
hampered by arbitrariness and uncertainties in the removal of the frictionally generated Joule heat. 
Building upon general pre-existing formulation, we implement a fully microscopic dissipation ap- 
proach which, based on a parameter-free, non-Markovian, stochastic dynamics, absorbs Joule heat 
cquivalently to a semi-infinite solid and harmonic substrate. As a test case, we investigate the stick- 
slip friction of a slider over a two-dimensional Lennard-Jones solid, comparing our virtually exact 
frictional results with approximate ones from commonly adopted dissipation schemes. Remarkably, 
the exact results can be closely reproduced by a standard Langevin dissipation scheme, once its 
parameters are determined according to a general and self-standing variational procedure. 



I. INTRODUCTION 

Ordinary, macroscopic sliding friction, a far reaching 
subject of enormous physical, technological and practi- 
cal importance, is notoriously complex and hard to ap- 
proach from a microscopic viewpoint, both experimen- 
tally and theoretically. The two last decades have seen 
quiet but important progress in that arena. Experimen- 
tally, the advent of nanosize slider methodologies is offer- 
ing much fresh data and lively progress. On the theory 
side, advances in computing hardware and codes now al- 
lows atomistic molecular dynamics (MD) simulations to 
be extensively used to describe sliding nanofriction: not 
simply as a mean of supplementing experimental studies, 
but as a general framework for gaining unique insight 
into the relevant tribologica l pr ocesses sometimes over- 
turning conventional wisdom^. In MD simulations, the 
classical dynamics of atoms is described by solving nu- 
merically Newton's equations of motion in a controlled 
computational experiment, where the interface geometry, 
sliding, boundary conditions and inter-particle interac- 
tions can be chosen and varied to explore various effects 
on friction, adhesion and wear. By following the parti- 
cle dynamics for a significant amount of time, quantities 
of physical interest such as instantaneous and average 
frictional force, mean velocities, heat flow, and correla- 
tion functions are calculated to characterize the sliding 
motion and the corresponding steady-state values. Un- 
like standard equilibrium MD simulations, friction mod- 
eling inherently involves dynamics and properties quite 
far from equilibrium. Moreover, as a rule, the dynamics 
is highly nonlinear too, for example in stick-slip friction. 

Actually, while MD simulations are quite valuable in 
qualitatively catching the physics of microscopic friction 
between extended solids, a quantitative agreement with 
experimental results is still beyond hopes ^. Besides 
the practical difficulty posed by the necessity to describe 
inter-atomic interactions by cither empirical force fields 
or with costly first principles calculations, an additional 



weak point of MD simulations lies in the impossibility to 
access the experimental time scaleiP. When attempting 
to simulate, e.g., a nanoscale Friction Force Microscopy 
experiment, with the tip advancing at a far low aver- 
age speed of ~ 1 /im/s, one can typically simulate a 
miserable ~ 1 pm advancement in a standard run, far 
too short to observe even a single atomic-scale event, let 
alone reaching a steady state, or the development of any 
instability process, and thus the quantitative evaluation 
of any useful frictional property. Therefore, whenever 
long-distance correlations and/or slow diffusive phenom- 
ena and/or long equilibration times are to be expected, 
fully atomistic MD approaches will only grab a qualita- 
tive scenario of the system tribological response. Nev- 
ertheless, there is so much direct physical insight to be 
extracted from MD simulations that it does make sense 
to run them even at larger speeds than in Atomic Force 
Microscopy (AFM) or Surface Force Apparatus experi- 
ments; and in fact, the sliding speed adopted in most cur- 
rent atomistic MD frictional simulations is much higher, 
in the 0.1 to 10 m/s range. 

The fast frictional motion in MD simulations ends up 
of course generating a vast amount of Joule heat. At 
the same time, the simulated system where that Joule 
heat is dispersed is generally of very limited size com- 
pared to the practically infinite environment of real ex- 
periments. That raises the problem, which is the focus 
of the present paper, of how that Joule energy can be 
continuously dissipated, "thermostated" away, in order 
for the simulated system to reach a realistic steady state 
rather than building up. At equilibrium, it does not mat- 
ter how the thermostat scheme is built, because equilib- 
rium properties do not depend on it. On the contrary, in 
dynamical non-equilibrium processes, such as those oc- 
curring in tribology under the action of external drive, 
the choice of a suitable physical thermostat is crucial, 
to dispose of the external energy which is continuously 
pumped into the system. In the framework of wearless 
friction, for instance, sliding-induced creation of phonons 
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is a crucial mechanism of energy dissipation. An unsolved 
problem in realistic MD is that the generated phonons 
cannot escape the small simulated contacting region be- 
tween a slider and the underneath substrate (see Fig. [T]) 
unlike in the real system, where they can properly dis- 
perse the Joule heat away from the interface. The simu- 
lation cell boundaries back-reflect the phonons towards, 
e.g., the slider-substrate contact, as shown in panel (b) 
of Fig. [2j affecting so the frictional response. As phonons 
are continuously generated by sliding, the simulated por- 
tions of the slider and substrate heat up, reaching quickly 
the melting point. Thus, in order to attain a frictional 
steady state in simulation, the Joule heat must be re- 
moved. Unfortunately, a realistic energy dissipation is 
generally impossible to mimic reliably, owing to size lim- 
itations of the simulation cell. The empirical introduction 
in the equations of motion of ad-hoc Langevin viscous 
damping terms —rrvyqi (with m and qi the mass and the 
velocity of the i-th substrate particle) and of an associ- 
ated random noise, corresponding to some "thermostat" 
temperature represents the handiest and commonest 
solution, which most simulations adopt. However, both 
this procedure and the choice of thermostat and damp- 
ing parameters 7 are vastly arbitrary. The problem is not 
just one of principle, for in many cases (including, just as 
a significant example, multiple-slips in AFIVP) the result- 
ing steady state and friction coefficient actually depend 
upon the choice of these unphysical parameters. Here, 
after demonstrating this unphysical dependence, we will 
pursue and detail a viable solution, whose core was al- 
ready outlined in a recent papei^. 




FIG. 1: Ideal block-scheme of a MD simulation of friction. 
To account properly for heat dissipation, the infinitely-thick 
substrate is divided into three regions: (i) a "live" slab com- 
prising layers whose atomic motion is fully simulated by New- 
ton's equations; (ii) a dissipative boundary zone, coincident 
with the deepmost simulated layer, whose dynamics includes 
effective damping (e.g., non-Markovian Langevin-type) terms, 
as in Eq. (pj; (iii) the remaining semi-infinite solid, acting as 
a heat bath, whose degrees of freedom are integrated out. 



II. NON-MARKOVIAN LANGEVIN 
APPROACH FOR REALISTIC TRIBOLOGICAL 
MODELING 

Basically, one wishes to modify the equations of motion 
inside a relatively small simulation cell so that they re- 
produce the frictional dynamics of a much larger system, 
once the remaining variables are integrated out. Inte- 
grating out degrees of freedom i s a tr aditional problem, 
largely analyzed in the literatur e! 5 * 8 * 9 !. In the context of 
MD simulation, Green's function methods were formu- 
lated for quasi-static mechanical contacts^; approaches 
based on a discrete-continuum matching have also been 
discussed^. Among others, time honored dissipation 
methods have been considered which replace the dynam- 
ics of the surrounding degrees of freedom (the "heat- 
bath" ) by several terms in the equations of motion for 
the system, describing effect gll^HMI sucn as 1) the renor- 
malization of the forces acting on and between the rele- 
vant coordinates; 2) the introduction of viscous drag de- 
scribing the energy dissipation from the system into the 
heat bath; 3) the introduction of random forces describ- 
ing the inverse effect of energy transfer from the bath 
into the system. Recently^ a direct implementation of 
a non-Markovian dissipation scheme, based on early for- 



mulations by Magalinskii and RubirJSLL 5 -! and s ubseq uent 
derivations by Li et alP^ and by KantoroviclP^i 8 ^, has 
demonstrated the correct disposal of friction-generated 
phonons in realistic MD simulations of sliding tribolog- 
ical systems, as the one sketched in Fig. [T] Once that 
was done, one could benchmark some simpler empirical 
Langevin scheme optimizing the 7 parameters so as to 
yield less arbitrary frictional properties. We describe 
here in detail how both goals are achieved, picking for 
our demonstration, without loss of generality, a specific 
two-dimensional (2D) realization. 

We consider a simplified tribological system red where 
the upper slider is represented by a one-dimensional 
(ID) chain of atoms along the x-axis driven on top 
of a 2D semi-infinite crystalline substrate lying in the 
(x, z) plane, where atoms interact, for simplicity, via 
first-neighbor Lennard- Jones (LJ) potential. The slider, 
pressed against the substrate by a normal "load" Fq, is 
driven along x (parallel to the surface) through a spring 
fc, whose end is pulled at constant velocity vq. Follow- 
ing earlier formulations^, the ideal infinitely thick sub- 
strate is divided, as sketched in in Fig. [T] in a 3D cartoon, 
into three regions: (i) an explicitly simulated substrate 
portion of N z atomic layers with displacement vectors 
r(t), (ii) the dissipative boundary layer, with displace- 
ment vectors q(t); and (iii) the remaining semi-infinite 
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FIG. 2: Propagation into the substrate of surface injected 
energy. Tapping on the surface layer (N z — 1) a burst of 
phonons has been created, its time evolution is monitored 
plotting the average kinetic energy of equi-spaced atomic lay- 
ers versus time, (a) shows a complete absorption of the 
phonon batch as it reaches the bottom of the simulation cell 
(Nz = 50) where our dissipation scheme is applied, (b) shows 
a total back-reflection of phonons when the correct dissipative 
kernels are switched off. 



solid acting as a phonon absorber, heat bath, with dis- 
placement vectors b(t). Under certain, not too restric- 
tive, assumptions described below, the heat bath degrees 
of freedom (iii) can be integrated out to let a small sim- 
ulation cell, namely (i)+(ii), account exactly for the en- 
ergy dissipation as due to a semi-infinite substrate, where 
the boundary layer (ii) is now ruled by effective non- 
Markovian Langevin equations, as derived in the follow- 
ing. The first needed assumption is to substitute the full 
LJ potential within regions (ii) and (iii), i.e. far away 
from the sliding interface, with its harmonic approxima- 
tion. This choice, necessary to derive an exact analyti- 
cal form for the effective forces acting on the boundary 
atoms, is all the more accurate the weaker the intensity 
of the slider perturbation and the lower the temperature. 
Nevertheless, for crystalline substrates well below the De- 
bye temperature, anharmonic perturbations reaching the 
heat bath can always be avoided by a sufficient thickness 
N z of the explicitly simulated substrate (i): these ex- 
citations, traveling through the LJ substrate, will grad- 



ually lose their energy turning into harmonic phonons 
prior approaching the boundary harmonic absorber. In a 
compact matrix notation, the hamiltonian of the system 
reads 

H = T + U{v, q) + qt . . q + qt . . b + • D • b, (1) 

where T is the overall kinetic energy term, U is the LJ in- 
teractions among atoms in region (i) and between regions 
(i) and (ii), 6 and d> are the LJ harmonic approximations 
for the atomic interactions in region (ii) and between re- 
gions (ii) and (iii) respectively, and D is the dynamical 
tensor of the heat bath (iii). Matrices and vectors have 
the form 
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where each component is again a matrix or a vector of 
components DJf„ or with latin indexes running over 
the atoms and greek indexes running over the two x and 
z coordinates. From the Hamiltonian ([T]), we can derive 
the following three sets of equations of motion: 
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dU(r,q) 
dx ' 
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-0 • q(t) - D • b(t) 
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(3) 

(4) 
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Notice that the dynamics of atoms in region (i) is influ- 
enced only by atoms of region (ii) , while the dynamics of 
region (iii) depends only upon the dynamics of region (ii) , 
in other words, thanks to the adopted cut-off LJ interac- 
tion, regions (i) and (iii) are decoupled and they interact 
only indirectly via the boundary layer (ii). The thick- 
ness size of the boundary layer (ii) depends on the cut-off 
radius: by considering only nearest-neighbors in the LJ 
interaction, we end up in our case with a region (ii) made 
of a single atomic layejf^l Thanks to the assumed har- 
monicity of the heat bath interactions, we can decouple 
the equations for b(t), diagonalizing the dynamical ten- 
sor D and finding its eigenvalues u>i and eigenvectors A;. 
By using the eigenvectors as a basis set b(i) = J2i 
we substitute this projection into Eq.([5]), obtaining a set 
of easily solvable decoupled equations for the normal co- 
ordinates £i, 
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where the rhs is a time dependent scalar quantity. The 
final expression for b(t) becomes 
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b(t) = A l • ( b (°) cos(cjit) + b(0) 
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0- / q(s) 



cos(cui(t — s)) 



ds A, 
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which depends on the initial conditions of atoms in region 
(iii) and on the actual position of atoms in region (ii). By 
substituting this expression into Eq.Q, we get 

mq(t) = - d ^ r ' q ) + (K(0)-g)-q(*)-m / K(t-s)q(s)ds- 
"Q Jo 

where K(t) and F(i) are defined as follows 
■(Aj-0)(0-Ai)" 



cos (ct^i), 



(9) 



f w = - • A *) A * ■ ( b (°) cos M) + b (°) 



sin(wii) 



(10) 



Equation ([8]) still depends on the initial conditions of the 
heat bath through F(t). Because this region is in princi- 
ple infinitely extended, we cannot specify the initial con- 
ditions for the position and the velocity of all its atoms; 
however, we are allowed to perform an equilibrium canon- 
ical ensemble average introducing a temperature T. Us- 
ing for the partition function the bath hamiltonian only, 
it is easy to prove that 
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being K B the Boltzmann's co nstant . Another possibil- 
ity, adopted for example in ref! 17 * 18 *, is to include in the 
partition function also the term ruling the interaction be- 
tween region (ii) and (iii). As a result the final effective 
equation of motion ( 16 1 takes a slightly different form. 



Using the previous conditions into Eq.(10), we end up 



with the following statistical properties for the force F(t) 

(F(t)) - 0, (F(t)F(f)) = mK B T K(t - t'), (14) 

or in component notation 

(F;{t)) = 0, (F;{t)Fl(t>)) = mk B TK%{t - f ). 

(15) 

Thus Eq.((8]) can be regarded as a non-Markovian 
Langevin equation with a gaussian random noise corre- 
lated according to the rules (15), and a dissipative term 
with a memory kernel function specified by ([9]). Its ex- 
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FIG. 3: Plot of some selected memory kernel functions versus 
time (LJ units). 



pression in single component notation is given by 
dU(v,q) 
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The first term takes into account the interaction between 
the boundary layer atoms and the rest of the simulated 
substrate. The second one is non-Markovian and non- 
conservative, introducing an effective damping propor- 
tional to the velocity of all the boundary layer atoms, 
via a time convolution with the me mory kernel func- 
tions K l ^ v if). The third term of Eq.(16) is the gaus- 



sian correlated noise ruled by the same memory ker- 
nel functions involved in the dissipation, in agreement 
with the fluctuation-dissipation theorem. Notice that in 
a standard Langevin equation the compliance with the 
fluctuation-dissipation theorem is imposed a priori and 
the noise properties are derived from this constraint. 
In our formulation, which starts from a microscopic set 
of Hamilton's equations, the fluctuation-dissipation the- 
orem is automatically fulfilled just performing the canon- 
ical ensemble average. The last term in Eq. ( 16 1 finally 
represents the harmonic coupling between i-th and j-th 
atoms within the boundary layer, where the coupling con- 
stant 9 l J u is modified by 1T^(0). This renormalization 
of the elastic coupling for the region (ii) atoms vanishes 
as we include the interaction between the bath and the 
boundary layer into the partion function for the ensemble 
average. It has been demonstrated theoretically-^, and 




FIG. 4: (a) Calculated friction force profile F(t) for the full non-Markovian dissipation scheme of Eq. \16\ , and for different 
empirical viscous damping schemes (b),(c) and (d) described in text, and identified by numbers 1-5 in Fig[5] Dashed lines: 
mean value (F) . 



it can be easily verified in simulations, that the applica- 
tion of Eq. ( [l6| to the boundary layer alone is sufficient 
to force the whole system to follow a canonical ensem- 
ble distribution with temperature T. The memory kernel 
matrix ([9]) in the single component notation reads 

r „ m x 



addition of more terms in the summations of (17 1. By 
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Each component is built from the harmonic eigenvalues 
and eigenvectors of the heat-bath dynamical matrix and 
from the coupling vectors <j) l ^ a containing the harmonic 
coupling constants of the i-th atom of region (ii) with 
the i-th heat-bath atom. As shown in Fig. |3j the ker- 
nels oscillate and decay rapidly with time, with power 
law tails due to the bath acoustical phonon branches. 
However as long as the heat bath region remains finite 
the summations in (17 1 are limited and the kernels are 
quasi-periodic functions 5 . Waiting for a large time A, 
which depends on the heat bath size, the kernel func- 
tions rise and decay again repeatedly, this time period- 
icity marks the energy back-reflections from one end of 
the finite heat bath to the opposite one. In the limit of 
infinitely extended heat bath A — > oo, no energy back 
reflection occurs. The numerical calculation of and 
Afc can be carried out only for a finite dynamical tensor, 
i.e. for a finite bath, however we can set Kjj v (t) = for 
all t > t with t < A preventing the first reflection. If 
the heat bath is large enough we verified that K l J u (t) for 
t < t is well converged, its shape being insensitive to the 



cutting kernels off after a time r one can limit the time- 
integrals in Eq. (16), which need to be evaluated at each 



time step, thus decreasing the heavy computational cost. 
But r represents also the maximum time for which the 
boundary layer retains memory and correlation, there- 
fore, via some convergence tests, we have to be sure that 
the quantities of interest do not depend on the chosen r 
value. Periodic boundary conditions along the x direc- 
tion guarantee translational invariance, so that K l J u (t) is 
a function of \i — j\ only. As kernels inherit their sym- 
metry properties from those of the heat-bath dynamical 
matrix, one can also show that K l 3(t) — K]? (t) and 
— Kv'JS)- When the separation \i — j\ grows, 
\K % i v {t)\ decrease, but again not exponentially, and cor- 
relations must be included up to large distance. Imple- 
menting this set of equations, along with ordinary New- 
ton's equations governing the remaining slider and sub- 
strate atom motion was our first MD simulation step. 
Figure[2ja) illustrates how a relatively thin (i.e. N z = 30 
layers) substrate (i+ii) is able to mimic the full ideal 
semi- infinite system (i+ii+iii) . Layer-resolved kinetic en- 
ergies inside the simulated substrate show a group of 
phonons initially created at the upper interface and prop- 
agating away from it. Upon reaching the boundary layer 
the phonons are perfectly absorbed as they propagate 
into the (integrated-out) semi- infinite crystal (iii). For 
comparison, Fig. [2jb) shows the same phonons massively 
back reflected once the memory kernels are removed from 
the boundary layer. 
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III. SIMULATING ATOMIC STICK-SLIP 

We next proceed to simulate sliding friction by driv- 
ing the slider (consisting, in the adopted 2D modeling, 
of a LJ chain of N' x = 9 atoms) over the live sub- 
strate, consisting of N z = 30 close packed layers and 
N x — 10 atoms per layer. Simulations were performed 
at temperature kgT = 0.035, roughly corresponding to 
T/T me iting = 0.06 (LJ units used throughout). To favor 
sliding, the strength of the slider-substrate LJ interaction 
is reduced from 1 to 0.6. The equations of motion are in- 
tegrated by a modified velocity- Verlet algorithm with a 
time step of At = 5 • 10~ 3 , and the memory kernel func- 
tions are cutoff at r = 5 ■ 10 3 time-stepsP3 Both the 
vertical load Fq and the lateral driving are applied to the 
slider center of mass, the equation of motion for the slider 
degrees of freedom Sj is 



dW(r,s) dU{s) 



d Six 

dU(r,s) 

dsi. 
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k{s xC M - v t) 



-Fq 



(18) 
(19) 



where U(r,s) is the LJ interaction with the substrate 
atoms (i) and U (s) is the L J interaction among the slider 
atoms, s x cm is the slider center of mass position along x. 
As usual the friction force is measured by the spring elon- 
gation k(s x CM ~ v ot) representing the slider resistence to 
the lateral driving. The applied load is Fo = 10, the av- 
erage sliding velocity vq = 0.01, and the spring constant 
k — 5. The result is the sawtooth force profile in Fig.Qa) 
typical of intermittent stick-slip friction. The friction 
coefficient, obtained by averaging over several stick-slip 
events, is (F)/F = 0.116 ± 0.002. The slider is slightly 
incommensurate with the substrate, so that the sawtooth 
pattern is quite irregular with a periodicity not exactly 
matching one lattice spacing. An anti-kink (physically 
corresponding to a tiny localized expansion in the parti- 
cle array density of the slider due to the interface mis- 
match) appears at the interface, moving in the opposite 
direction with respect to the slider: the height of the 
sawtooth spikes is proportional to the jump length of the 
anti-kink. Higher spikes occur for simultaneous forward 
jumps of many atoms, smaller ones correspond to jumps 
of 2 — 3 atoms at once. A measure of the distribution of 
the spike heights is the variance of F(t), i.e., 



1 



[F(t) - (F)) 2 dt, 



(20) 



where t s is the total simulation time. Numerical simu- 
lations carried out with the full Eq. (16 1, and the cor- 



responding frictional results are essentially exact for the 
system considered. That completes our first important 
goal of implementing the correct Joule heat removal, thus 
also establishing a benchmark reference. Not surpris- 
ingly, this numerical implementation is time consuming. 
In particular the computational effort required to inte- 
grate the non-Markovian term, where boundary atoms 
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FIG. 5: (a) and (b) illustrate the friction coefficient (F)/Fq 
and variance (a) behaviors as a function of the damping coeffi- 
cient 7 for different empirical Langevin dissipation schemes, in 
comparison with the exact values from the full non-Markovian 
simulation (gray stripes), (c) shows the boundary layer ab- 
sorbed energy W of the Langevin thermostat (211. Note the 



good coincidence of exact and empirical frictional behavior 
for the optimal 7 that maximizes W. 



are strongly correlated, scales as N x . Carrying out fu- 
ture fully realistic frictional simulations for large-size 3D 
sliding systems within this scheme is in our view entirely 
possible, but may pose some practical challenge of par- 
allel computing. 

This brings us to our second point. As was men- 
tioned, much simpler and faster approximate frictional 
simulations are realized once the non-Markovian mem- 



ory kernels of Eq. ( 16 ) are empirically replaced with 



a more ordinary Markovian Langevin viscous damping 
— rwyql (t), along with the appropriate gaussian stochas- 
tic force Ri(t) with (B^(t)) = and (i^(t)i%(i')) = 
2mkBT r )8 u >u 5i jS(t — t'), so that the equation of motion 
of the z-th thermostated atom in the system reads 



mqUt) 



dU 



(21) 



where Li is the LJ inter-atomic interaction. Performing 
a series of simulations with the same system parameters, 
the previous exact implementation now offers the possi- 
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bility to benchmark the empirical damping 7. In prin- 
ciple, this standard Langevin scheme can be differently 
exploited, applying it to: I. (Figj4fb), curves 1 and 2) the 
slider atoms only, while freezing the substrate degrees of 
freedom, as typically done in the simplified framework of 
Prandtl-Tomlinson and Frenkel-Kontorova modeling 19 ; 
II. (FigQc), curves 3 and 4) to each substrate atom, pos- 
sibly by making it site-dependentP^ III. (FigQd), curve 
5) just to the bottom simulation-cell layer, as consid- 
ered for the parameter-free, non-Markovian, stochastic 
dynamics. In all these cases, we find a strong depen- 
dence of the system frictional response on the choice of 
the damping 7, in general deviating always systemati- 
cally from the correct benchmark. Figure ^ a), (b) shows 
the behavior of the friction coefficient and its variance, 
respectively, as a function of 7. The grey stripes indi- 
cates the benchmark values of (F ) / Fo and (a) obtained 
with our parameter-free dissipation scheme, mimicking a 
semi-infinite substrate. The dashed line represents the 
results for standard Langevin equations applied, only, to 
the slider atoms (case I.): this turns out to be the most 
unrealistic and 7-sensitive situation. A too large 7 in- 
troduces a strong viscous character, and leads to over- 
estimating the friction force, while a too small 7 results 
in a chaotic behavior, with the slider dynamics being un- 
able to dissipate enough energy. This scenario is outlined 
in Fig. [l|b) where the F(t) stick-slip profile is plotted 
for 7 = 0.01 (blue solid line) and for 7 = 1.0 (black 
dotted line). At 7 = 0.035 in Fig. [5](a), this average 
friction force curve crosses the "exact-method" (grey) 
stripe with a value (F)/F = 0.117 ± 0.001, but with 
a too large variance (a = 0.28) (Fig. [5jb) ) , and a conse- 
quent very inaccurate reproduction of the stick-slip pat- 
tern (not shown). The dotted line in Fig. [5] represents 
the Markovian Langevin thermostat applied, more real- 
istically, to all substrate atoms (case II.): the slider ex- 
changes energy with the substrate by exciting phonons at 
the interface; these phonons are then damped within the 
substrate independently of the slider velocity. However, 
a too large 7 will lead to a very viscous surface prevent- 
ing the correct energy exchange between the slider and 
the substrate, and (F) / Fo increases too much, as in the 
previous case. A too small 7, on the contrary, makes 
the substrate unable to dissipate the phonons, which are 
then reflected back, reaching again the surface and heat- 
ing it to unphysically large temperatures, thus spuriously 
decreasing the friction force. Fig. [4jc) , corresponding to 
such case II., shows F(t) for a low 7 value of 0.01 (blue 
solid line): the effect of the reflected phonons is to reduce 
the static friction force, decreasing the swing of the saw- 
tooth profile. F(t) is also displayed for 7 = 0.1 (black 
dotted line): the average friction force here approaches 
our semi- infinite substrate result, mimicking well also the 
stick-slip profile, as highlighted by the simultaneous good 
values of the friction coefficient and the standard devia- 
tion in Fig.|5ja),(b). However, there is here (case II.) no a 
priori possibility to choose the optimal value of the damp- 
ing parameter without having previously performed an 



exact non-Markovian benchmark calculation. Besides, in 
order not to directly interfere with the detailed dynamics 
and the slider-substrate energy exchange, the Langevin 
viscous damping term should be switched on far from 
the surface as, e.g., in the bottom dissipation layer (case 
III.), shown by the continuous line in Fig. [5] We find that 
there exists an optimal damping 7 opt (here 7 opt ~ 10) for 
which both the friction coefficient and its variance agree 
well with the exact values (see Fig. [5^ a), (b)). Moreover, 
also the stick-slip profile in Fig. [I]for 7 = j opt (panel(d)) 
compares excellently with the exact one (panel(a)). Re- 
markably the 7 value for which the friction profile bet- 
ter resemble the exact one corresponds to the one which 
maximizes the average friction force. In order to under- 
stand this relation, we look at the energy dissipated by 
the boundary layer: 

W = -m^ J 7 q l • def = -m^2 J 7 W?dt, (22) 

i i 

finding a maximum at the same 7 values as illustrated in 
Fig.gc). This maximum occurs because back-reflection 
of phonons is large both when the boundary layer damp- 
ing 7 is too small and too large. The efficiency in the 
energy removal goes as —rrvyqi, so that at low 7 values 
the boundary layer atoms cannot dissipate significantly 
even if vibrating very fast; in the opposite limit of large 
7, the boundary layer dynamics becomes so viscous (low 
atomic velocities) that an effective dissipation is again 
hampered. At 7 = 7 op t, we reach a good compromise be- 
tween the strength of the damping and the atom veloci- 
ties and most of the impinging energy is disposed of. The 
agreement between the exact frictional results, where no 
phonons are back reflected, and the approximate ones 
is therefore best when energy back reflection is minimal 
and this can occur for a single 7 value only. The mini- 
mal phonon back reflection condition also establishes the 
smallest temperature at the sliding interface. While this 
makes good physical sense, we still contemplate the possi- 
bility that the numerical result might be just some kind of 
coincidence in a single simulation. We therefore proceed 
to change system parameters, including sliding velocity, 
and load. In all cases we find an optimal 7 value, where 
both the friction force and the energy dissipated by the 
boundary layer are maximized and where both average 
friction and variance coincide with the exact value sepa- 
rately calculated by a full non-Markovian simulation. For 
example the variable load results of Fig. [6] show that the 
coincidence of optimal and exact friction is systematic 
as well as the presence of the force maximum that can 
be thus exploited as a tool to calibrate the viscous coef- 
ficient 7 for any general system even without the exact 
non-Markovian benchmark. 



IV. CONCLUSIONS 

We have shown here that sliding friction obtained by 
Molecular Dynamics simulations may depend heavily on 
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FIG. 6: Average friction force for different loads Fo. Gray 
stripes show the values obtained with the non-Markovian ap- 
proach in comparison with the black curves for the boundary 
Langevin scheme at different 7. 



by the stick-slip profiles in Fig. [2j More importantly, 
these schemes generally offer no clue on how to optimize 
the empirical parameter 7 in the absence of the exact 
simulation. 

We then showed how the real dissipation of phonons 
into a harmonic semi-infinite solid substrate can be sim- 
ulated by implementing well established non-Markovian 
schemes. Once the exact non-Markovian dissipation is re- 
placed by an approximate and empirical Langevin damp- 
ing 7 applied to the bottom layer of the simulated sub- 
strate slab, an optimal value for 7 is easily and varia- 
tionally found by maximizing dissipation - a condition 
which can be established without resort to any exact ref- 
erence calculation. This is a result which in all likeli- 
hood appears more general than the simple model used 
to demonstrate it, and should thus be quite valuable for 
general applications. 



the scheme adopted for the elimination of Joule heat. 
None of the empirical but commonly used dissipation 
schemes seems satisfactory. One might for example ap- 
ply a Langevin viscous damping 7 to the slider atoms 
alone^l, or, uniformly to all substrate atoms^. Shown 
as dashed and dotted lines respectively in Fig. [3j the 
friction coefficients produced by these approximations, 
although crossing the correct values as a function of 7, 
generally yield a much lower quality description as seen 
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